Preparación del entorno de trabajo y carga de librerías
# Limpieza
rm(list = ls())
suppressPackageStartupMessages({
library(data.table)
library(dplyr)
library(caret)
library(scales)
library(ggplot2)
library(stringi)
library(stringr)
library(dataPreparation)
library(knitr)
library(kableExtra)
library(ggpubr)
library(tictoc)
library(ggeasy)
library(lubridate)
library(inspectdf)
library(ranger)
library(gbm)
})Utilizando los datos provistos por el Ministerio del agua de Tanzania, se requiere construir un modelo que sea capaz de predecir cuales bombas de agua están operativas, operativas pero que necesitan reparación o están dañadas, basadas en un set de datos train.
De forma original se incluye los dataset train y test mas el archivo con las labels.
# Data Loading
#Fichero datos train
vtrain <- fread("train_set.csv")
vtrain$flag <- 1 # Columna que indica si es parte del set train (1) o test (0)
#Fichero datos test
vtest <- fread("test_set.csv")
vtest$flag <- 0 # Columna que indica si es parte del set train (1) o test (0)
#Fichero con labels o objetivo
vlabels <-fread("labels.csv")
# Se unen las labels con el set de datos train
train <- merge(vlabels, vtrain)
# Se unen ambos datasets (train y test) lo cual es lo recomendado para mas adelante trabajar en Feature engineering
datos <- as.data.table(rbind(vtrain, vtest))
#Comprobacion
head(datos)head(vlabels)#Distribucion de los datos
str(datos)## Classes 'data.table' and 'data.frame': 74250 obs. of 41 variables:
## $ id : int 69572 8776 34310 67743 19728 9944 19816 54551 53934 46144 ...
## $ amount_tsh : num 6000 0 25 0 0 20 0 0 0 0 ...
## $ date_recorded : IDate, format: "2011-03-14" "2013-03-06" ...
## $ funder : chr "Roman" "Grumeti" "Lottery Club" "Unicef" ...
## $ gps_height : int 1390 1399 686 263 0 0 0 0 0 0 ...
## $ installer : chr "Roman" "GRUMETI" "World vision" "UNICEF" ...
## $ longitude : num 34.9 34.7 37.5 38.5 31.1 ...
## $ latitude : num -9.86 -2.15 -3.82 -11.16 -1.83 ...
## $ wpt_name : chr "none" "Zahanati" "Kwa Mahundi" "Zahanati Ya Nanyumbu" ...
## $ num_private : int 0 0 0 0 0 0 0 0 0 0 ...
## $ basin : chr "Lake Nyasa" "Lake Victoria" "Pangani" "Ruvuma / Southern Coast" ...
## $ subvillage : chr "Mnyusi B" "Nyamara" "Majengo" "Mahakamani" ...
## $ region : chr "Iringa" "Mara" "Manyara" "Mtwara" ...
## $ region_code : int 11 20 21 90 18 4 17 17 14 18 ...
## $ district_code : int 5 2 4 63 1 8 3 3 6 1 ...
## $ lga : chr "Ludewa" "Serengeti" "Simanjiro" "Nanyumbu" ...
## $ ward : chr "Mundindi" "Natta" "Ngorika" "Nanyumbu" ...
## $ population : int 109 280 250 58 0 1 0 0 0 0 ...
## $ public_meeting : logi TRUE NA TRUE TRUE TRUE TRUE ...
## $ recorded_by : chr "GeoData Consultants Ltd" "GeoData Consultants Ltd" "GeoData Consultants Ltd" "GeoData Consultants Ltd" ...
## $ scheme_management : chr "VWC" "Other" "VWC" "VWC" ...
## $ scheme_name : chr "Roman" "" "Nyumba ya mungu pipe scheme" "" ...
## $ permit : logi FALSE TRUE TRUE TRUE TRUE TRUE ...
## $ construction_year : int 1999 2010 2009 1986 0 2009 0 0 0 0 ...
## $ extraction_type : chr "gravity" "gravity" "gravity" "submersible" ...
## $ extraction_type_group: chr "gravity" "gravity" "gravity" "submersible" ...
## $ extraction_type_class: chr "gravity" "gravity" "gravity" "submersible" ...
## $ management : chr "vwc" "wug" "vwc" "vwc" ...
## $ management_group : chr "user-group" "user-group" "user-group" "user-group" ...
## $ payment : chr "pay annually" "never pay" "pay per bucket" "never pay" ...
## $ payment_type : chr "annually" "never pay" "per bucket" "never pay" ...
## $ water_quality : chr "soft" "soft" "soft" "soft" ...
## $ quality_group : chr "good" "good" "good" "good" ...
## $ quantity : chr "enough" "insufficient" "enough" "dry" ...
## $ quantity_group : chr "enough" "insufficient" "enough" "dry" ...
## $ source : chr "spring" "rainwater harvesting" "dam" "machine dbh" ...
## $ source_type : chr "spring" "rainwater harvesting" "dam" "borehole" ...
## $ source_class : chr "groundwater" "surface" "surface" "groundwater" ...
## $ waterpoint_type : chr "communal standpipe" "communal standpipe" "communal standpipe multiple" "communal standpipe multiple" ...
## $ waterpoint_type_group: chr "communal standpipe" "communal standpipe" "communal standpipe" "communal standpipe" ...
## $ flag : num 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, ".internal.selfref")=<externalptr>
Se despliega la distribución de la variable objetivo (contenida en labels) Se observa algo de desbalance el cual no debería afectar demasiado
cuenta <- vlabels %>% count(status_group)
porcentaje <- round( prop.table(table(vlabels$status_group))*100, 2)
kable(cuenta, col.names = c('status_group', 'count'))| status_group | count |
|---|---|
| functional | 32259 |
| functional needs repair | 4317 |
| non functional | 22824 |
kable(porcentaje, col.names = c('status_group', '%'))| status_group | % |
|---|---|
| functional | 54.31 |
| functional needs repair | 7.27 |
| non functional | 38.42 |
barplot(porcentaje, col=rgb(0.2,0.4,0.6,0.6))# categorical plot
x <- inspect_cat(datos)
show_plot(x)# correlations in numeric columns
x <- inspect_cor(datos)## Warning: Columns with 0 variance found: flag
show_plot(x)# feature imbalance bar plot
x <- inspect_imb(datos)
show_plot(x)# memory usage barplot
x <- inspect_mem(datos)
show_plot(x)# missingness barplot
x <- inspect_na(datos)
show_plot(x)# histograms for numeric columns
x <- inspect_num(datos)
show_plot(x)# barplot of column types
x <- inspect_types(datos)
show_plot(x)Una gran parte de las features son de tipo categóricas
Se deben explorar los missing values (en gris) en el primer gráfico de frecuencias de las categóricas.
wpt_name, subvillage, scheme_name e installer son las que ocupan mas espacio en memoria, si se revisan en el gráfico de frecuencias de las categóricas, son las que mas categorías contienen por cada una de las variables.
public_meeting y permit muestran columnas con % de NA (5.6% y 5.1% respectivamente)
La feature construction_year tiene valores en 0, es decir, no hay información codificada del año de construcción de la bomba de agua.
Se construye un dataset sencillo para comenzar, que contiene solo variables numéricas. A medida que se avanza en el conocimiento del set de datos se irá agregando complejidad.
datos_num <- Filter(is.numeric, datos)
head(datos_num)Observando las variables numéricas hay algunas que pueden resultar interesantes y otras que deberían ser descartadas del estudio.
Para comprobar lo anterior, en primer lugar se estudiara la presencia de valores en 0 que probablemente correspondan a valores NA.
# Syntax
temp <- datos_num %>% mutate_at(vars(-c(flag)), ~na_if(., 0))
x <- inspect_na(temp)
show_plot(x)Basado en lo anterior y dado el alto num de NA las variables num_private y amount_tsh saldrán del modelo.
temp <- temp %>% select(-num_private, -amount_tsh)Se convierten los na anteriores en ceros otra vez.
temp[is.na(temp)] <- 0
head(temp)Se construyen algunas variables nuevas que pueden ser útiles a la hora de modelar. Se agregan con una fe_ para identificarlas.
La variable construction_year tiene un % importante de valores missing. Se crea una columna nueva con valores missing imputados segun la media.
temp$fe_construction_year<-round(ifelse(temp$construction_year==0, mean(temp$construction_year[temp$construction_year>0]),temp$construction_year), 0)#Combinacion de latitud y longitud
temp$fe_lonlat <- sqrt(temp$longitude^2 + temp$latitude^2)Otra variable que calcula la antigüedad de la bomba basada en su fecha de construcción
year <- year(now())
temp$fe_antiguedad <- (year - temp$fe_construction_year)Se elimina la variable construction_year
temp$construction_year <- NULL# Separa train y test segun su flag
trainset<-temp[temp$flag==1,]
#table(trainset$flag) comprobacion
testset<-temp[temp$flag==0,]
# Se combina el train set con la variable objetivo
trainset <- merge(trainset, vlabels, by ='id', sort = FALSE)
# Elimina la columna flag y la columna id
trainset$flag <- NULL
testset$flag <- NULL
trainset$id <-NULLLa columna status_group indica en palabras si es funcional o no. Se modifica para hacer mas simple la lectura
trainset = trainset %>%
mutate(status_group = ifelse(status_group== "functional", 0,
ifelse(status_group == "functional needs repair",1 , 2)))
table(trainset$status_group)##
## 0 1 2
## 32259 4317 22824
Tomando como base lo anterior y solo considerando variables numéricas se construye el primer modelo
fit <- ranger(status_group ~. ,
data = trainset,
num.trees = 300,
mtry=6,
importance = 'impurity',
write.forest = TRUE,
min.node.size = 1,
splitrule = "gini",
verbose = TRUE,
classification = TRUE,
seed=1234
)Se despliegan los resultados
print(fit)## Ranger result
##
## Call:
## ranger(status_group ~ ., data = trainset, num.trees = 300, mtry = 6, importance = "impurity", write.forest = TRUE, min.node.size = 1, splitrule = "gini", verbose = TRUE, classification = TRUE, seed = 1234)
##
## Type: Classification
## Number of trees: 300
## Sample size: 59400
## Number of independent variables: 9
## Mtry: 6
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: gini
## OOB prediction error: 28.35 %
Se crea la matriz de confusión para trainset
predictions_train <- predict(fit, data = trainset)
confusionMatrix(table( trainset$status_group, predictions_train$predictions))## Confusion Matrix and Statistics
##
##
## 0 1 2
## 0 32189 0 70
## 1 362 3931 24
## 2 386 0 22438
##
## Overall Statistics
##
## Accuracy : 0.9858
## 95% CI : (0.9848, 0.9868)
## No Information Rate : 0.5545
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9741
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Statistics by Class:
##
## Class: 0 Class: 1 Class: 2
## Sensitivity 0.9773 1.00000 0.9958
## Specificity 0.9974 0.99304 0.9895
## Pos Pred Value 0.9978 0.91059 0.9831
## Neg Pred Value 0.9724 1.00000 0.9974
## Prevalence 0.5545 0.06618 0.3793
## Detection Rate 0.5419 0.06618 0.3777
## Detection Prevalence 0.5431 0.07268 0.3842
## Balanced Accuracy 0.9873 0.99652 0.9927
Se predice sobre los datos del concurso
predictions_concurso <- predict(fit, data = testset)
resultados_concurso<- as.data.frame(cbind( testset$id,predictions_concurso$predictions))
names(resultados_concurso)<-c("id", "status_group")
resultados_concurso$status_group<-ifelse(resultados_concurso$status_group==0,"functional", ifelse(resultados_concurso$status_group==1,"functional needs repair","non functional"))
# Crea archivo para subir a DataDriven
fwrite(resultados_concurso, file = "results_model1.csv")Variables importantes
vars_imp <- fit$variable.importance
vars_imp <- as.data.frame(vars_imp)
vars_imp$myvar <- rownames(vars_imp)
vars_imp <- as.data.table(vars_imp)
setorder(vars_imp, -vars_imp)Plot de variables más importantes
ggbarplot(vars_imp,
x = "myvar", y = "vars_imp",
#fill = 'myvar',
color = "blue",
palette = "jco",
sort.val = "asc",
sort.by.groups = FALSE,
x.text.angle = 90,
ylab = "Importancia",
xlab = 'Variable',
#legend.title = "MPG Group",
rotate = TRUE,
ggtheme = theme_minimal()
)Las features con mas peso predictivo para este modelo son latitud, la transformación de longitud y latitud, longitud y gps_height
Este primer modelo tiene un scoring de 0.7136 lo cual no es el más optimo pero para ser un modelo de tipo básico (como punto de partida y solo con datos numéricos) está aceptable. El próximo modelo incluirá variables numéricas como categóricas y algunas transformaciones para esas variables, buscando mejorar el scoring.